Leaflet Homework Key

library(dplyr)
library(sf)
library(tidyverse)
library(leaflet)
library(sf)
library(leaflet.providers)
library(RColorBrewer)
library(tigris)

Dot Density Plot

education_url <- "https://github.com/dyerlab/ENVS-Lectures/raw/master/data/edu_attainment-shp.zip"

file <- tempfile(pattern = "edu", tmpdir = ".", fileext = ".zip")
file

[1] ".\\edu253c7a4044d8.zip"

download.file(education_url, destfile = file)

trying URL 'https://github.com/dyerlab/ENVS-Lectures/raw/master/data/edu_attainment-shp.zip'
Content type 'application/zip' length 285435 bytes (278 KB)
downloaded 278 KB

unzip(file, exdir = ".")


points <- st_read("edu_attainment-shp")

Reading layer `edu_attainment' from data source `C:\Users\Charis\Documents\Spring2021\GISinR\GISinR\ENVS-Lectures\edu_attainment-shp' using driver `ESRI Shapefile'
Simple feature collection with 15745 features and 1 field
geometry type:  POINT
dimension:      XY
bbox:           xmin: -77.59891 ymin: 37.44779 xmax: -77.38568 ymax: 37.60228
geographic CRS: WGS 84

Convert education levels to factor - this allows us to use colorFactor() to create our color palette

points$education <- factor(points$education, levels = c("No HS diploma", "HS, no Bachelors", "Bachelors" , "Post-Bachelors"  )) 

Create color palette

factpal <- colorFactor(brewer.pal(n = 4, name = "RdBu"), points$education) #used to color categorical data

Create dot density plot

leaflet(points) %>% 
  addProviderTiles(providers$CartoDB.Positron)%>% #add provider tiles
  addCircles(stroke = FALSE,  fillOpacity = 1, #add points
              color = ~factpal(education)) %>% #color based on educational attainment
  addLegend(pal = factpal, values = ~education, title = "Education Attainment")%>% #add legend
  addScaleBar(position = "bottomright") #add scale bar

Choropleth plot

load in the census tract data

tracts <- tracts("VA", "Richmond city") %>%
  st_transform(4326)

Using FIPS code '51' for state 'VA'
Using FIPS code '760' for 'Richmond city'

Data Prep

edu_tracts <- tracts %>%
  st_buffer(dist = 0) %>% # buffer by 0 to avoid error message
  st_intersection(points) # intersect the two data sets

edu_tracts %>%
  group_by(NAME, education) %>% # group by tract name and education level
  tally() -> tally # tally() counts the number of rows for the groups, so output is one row for every attainment level for each tract (so if a tract has all 4 attainment levels, there are 4 rows of data for that tract)

tally

Simple feature collection with 246 features and 3 fields
geometry type:  GEOMETRY
dimension:      XY
bbox:           xmin: -77.59891 ymin: 37.44779 xmax: -77.38568 ymax: 37.60228
geographic CRS: WGS 84

tally <- as.data.frame(tally) # converted to data.frame to be able to drop the geometry because the `pivot_wider()` function would not work the way I wanted it to when the geometry was present

tally <- tally %>%
  select(NAME, education, n)%>% # do not want the geometry 
  pivot_wider(names_from = education, values_from = n) # pivot the table so that instead of having multiple rows per tract, you have one row per tract with a column for each attainment level

tally[is.na(tally)] <- 0 # define NA values as zero

Create attainment metric for choropleth

tally <- tally %>%
  mutate(edu = (((`No HS diploma` * 1) + (`HS, no Bachelors` * 2) + (`Bachelors` * 3) + (`Post-Bachelors` * 4))/
                  (`No HS diploma`  + `HS, no Bachelors` + `Bachelors` + `Post-Bachelors`)))

# I created a column, edu, and did a weighted average for each row. So 'No HS diploma' was weighted as 1, 'HS, no Bachelors' as 2, ect. 

# I used this new column as my average education attainment value to color my polygons

tally %>%
  select(NAME, edu) -> tally # no longer need the individual count data

Join the education data back with the tract data

tracts %>%
  left_join(tally, by = "NAME") -> tracts

Create the palette

pal <- colorNumeric(
  palette = "RdBu",
  domain = tracts$edu) # This time I want a continuos palatte, so I used `colorNumeric()`

Make the choropleth map

leaflet(tracts) %>%
  addProviderTiles(providers$CartoDB.Positron)%>% # add provider tiles
  addPolygons( fillOpacity = 1, # add tract polygons
               fillColor = ~pal(edu), smoothFactor = 0.2, # fill color based on the calculated 'edu' metric
               color = "grey", 
               opacity = 1,
               weight = 1) %>%
  addLegend(pal = pal, values = ~edu, title = "Education Attainment") # add legend

I did not like this legend, as it’s unclear what the numbers mean. With some googling I found a solution

labels <-c("No HS diploma (1)", "HS, no Bachelors (2)",
           "Bachelors (3)" , "Post-Bachelors (4)"  ) # define the labels for the legend

leaflet(tracts) %>%
  addProviderTiles(providers$CartoDB.Positron)%>% # add provider tiles
  addPolygons( fillOpacity = 1, # add tract polygons
               fillColor = ~pal(edu), smoothFactor = 0.2, # fill color based on the calculated 'edu' metric
               color = "grey", 
               opacity = 1,
               weight = 1) %>%
  addLegend(pal = pal, values = ~edu, 
            opacity = 1,
            title = "Education Attainment",
            labFormat = function(type, cuts, p) {  # Here's the trick
              paste0(labels)
            })

That looks better, I think

Combine the maps

Final map